A numerical testbed for singularity excision in moving black hole spacetimes 
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We evolve a scalar field in a fixed Kerr-Schild background geometry to test simple (3+l)-dimensional 
algorithms for singularity excision. We compare both centered and upwind schemes for handling 
the shift (advection) terms, as well as different approaches for implementing the excision boundary 
conditions, for both static and boosted black holes. By first determining the scalar field evolution in 
a static frame with a (1 + 1) -dimensional code, we obtain the solution to very high precision. This 
solution then provides a useful testbed for simulations in full (3+1) dimensions. We show that some 
algorithms which are stable for non-boosted black holes become unstable when the boost velocity 
becomes high. 



I. INTRODUCTION 

The long-term numerical evolution of black holes is 
one of the most important and challenging problems 
in numerical relativity. Simultaneously, it is a prob- 
lem for which a solution is very urgently needed; binary 
black holes are among the most promising sources for 
the gravitational wave laser interferometers currently un- 
der development, including LIGO, VIRGO, GEO, TAMA 
and LISA, and theoretically predicted gravitational wave 
templates are crucial for the identification and interpre- 
tation of possible signals. 

Numerical difficulties arise from the complexity of Ein- 
stein's equations and the existence of a singularity in- 
side the black hole (BH). Numerical simulations based on 
the traditional Arnowitt-Deser-Misner (ADM) decompo- 
sition in (3 + 1) dimensions, for example, often develop 
instabilities . The gauge (coordinate) freedom inher- 
ent to general relativity constitutes a further complica- 
tion. Singularity avoiding slicings (^-^) can follow evolu- 
tions involving black holes only for a limited time, since 
the stretching of time slices typically causes simulations 
to crash on time scales far shorter than the time required 
for a binary BH orbital period. 

Recently, there have been several very promising nu- 
merical breakthroughs. Stable formulations of Ein- 
stein's equations using a conformal-tracefree decompo- 
sition have been developed by Shibata and Nakamura 
H and by Baumgarte and Shapiro The so-called 
BSSN formulation has shown remarkable stability prop- 
erties when compared to the ADM formulation in a wide 
range of numerical simulations [0-^2|. Also, the develop- 
ment of "singularity excision" techniques Jl3[-p7t, excis- 
ing the singularity from the computational domain, may 
allow for long term binary BH evolutions (see |l8| ] for pre- 
liminary results). There have also been proposals 19 2l]] 
for alternative families of binary initial data based on the 
Kerr-Schild form of the Schwarzschild (or Kerr) metric to 
represent each of the BHs. The use of Kerr-Schild coor- 
dinates is desirable for the numerical evolution of BHs 



and suitable for applying singularity excision techniques 
because the coordinates smoothly penetrate the horizon 
of the holes. Moreover, they provide a natural frame- 
work for constructing initial data without assuming con- 
formal flatness, essential for representing non-distorted 
Kerr Black holes. 

Typically, the numerical grid extends into the black 
hole in singularity excision applications. Moreover, the 
shift vector sometimes becomes larger than unity, so 
that some finite difference stencils become "acausal" and 
therefore unstable. To avoid this problem, "causal dif- 
ferencing" [[l4||l(f | or "causal reconncction" schemes 0] 
have been suggested (see also @,|3 - 28 ) . Unfortunately, 
these schemes are fairly involved and complicated. 

Alcubierre et al p9| , |30[ | recently proposed a simple BH 
excision algorithm for a non-boosted distorted BH evo- 
lution which avoids the complications of causal differenc- 
ing. Their algorithm is based on the following simpli- 
fications: (a) excise a region adapted to Cartesian co- 
ordinates, (b) use a simple inner boundary condition at 
the boundary of the excision zone, (c) use centered dif- 
ferences in all terms except the advection terms on the 
shift, where upwind differencing along the shift direction 
is used. Their algorithm is quite successful and allows for 
accurate and stable evolution of non-boosted distorted 
BHs for hundreds of dynamical times. 

Here we devise an experiment to test some simple sin- 
gularity excision algorithms for evolving dynamical fields 
in numerical 3D black hole spacetimes. Our goal is to 
identify a numerical scheme which is simple and stable 
like the one presented in |2^] for non-boosted BHs but 
can also handle moving BHs. Our experiment consists 
of tracking the propagation of a scalar field in the fixed 
background spacetime of a Kerr-Schild BH. By first solv- 
ing for the evolution in a static frame by means of a ID 
code, we obtain the solution to arbitrary accuracy. We 
then perform simulations in full 3D for both stationary 
and boosted Kerr-Schild BHs, using our ID "exact" re- 
sults for detailed comparison. We use this testbed to 
compare schemes which handle advection terms by cen- 
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tered versus upwind differencing. We also run our experi- 
ment for different implementations of the excision bound- 
ary conditions. Our aim is to find numerical stencils 
which we can later adapt to a BSSN scheme for evolving 
dynamical spacetimes with moving (e.g., binary) BHs. 
While the proving ground we construct here employs a 
fixed background geometry, we feel that it furnishes nec- 
essary, if not sufficient, conditions that must be met for 
a field evolution scheme in any numerical spacetime con- 
taining BHs. We provide this report to convey the utility 
of this testbed as a quick diagnostic of alternate differ- 
encing stencils and with the hope that it might prove 
helpful for other 3D code builders. 

Our results show that, in general, an upwind scheme 
is more stable than a centered scheme. This is consis- 
tent with the results of p9[| . We also find that a higher 
resolution is needed for an upwind than for a centered 
scheme to achieve a desired accuracy, due to the diffu- 
sive character of an upwind scheme. We find that the 
numerical implementation used in ]29| ] is stable in the 
non-boosted case but unstable in the boosted case. How- 
ever, the stability can be restored by using an alternative 
excision boundary condition which is stable in both the 
non-boosted and boosted cases. In all the cases studied, 
the use of a third-order extrapolation condition at the 
excision boundary is required for stable runs. 

The paper is organized as follows: We describe the 
Kerr-Schild BH spacetime and the scalar field equation 
in this background geometry in Sec. |fi|. Sec. Ill is devoted 
to a discussion of initial data, numerical algorithms, and 
different boundary conditions. We present our 3D nu- 
merical results, and compare them to very accurate ID 
results in Sec. IV. We summarize and discuss the im- 
plications of our findings in Sec. fv| We also include two 
Appendices. Appendix [A| sketches the von Neumann sta- 
bility analysis of the ID centered and upwind schemes. 
Appendix [Bj describes the Lorentz transformation of a 
scalar wave. Throughout the paper we adopt geometrized 
units with G = c = 1 . 



t he horizon are well-behaved. Comparing the metric 
(pT| ) with the "ADM" metric typically used in 3 + 1 for- 
mulations, one identifies the lapse function a, shift vector 
fti and the spatial 3-metric ^'g%j as 



A = 2Hl t £i, 
(3) 5y = Vv+ZHtd 



2H£l 



(2.3) 



For the time-independent Schwarzschild spacetime (the 
"Eddington-Finkelstein form" |32|]) in Cartesian coordi- 
nates we have 
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where M is the total mass-energy and r 2 = (x 1 ) 



(x 2 ) 2 + (x 3 ) 2 . The Kerr-Schild metric lyj) is form- 
invariant under Lorentz transformations. Applying a 
constant Lorentz transformation A (with boost velocity 
v a s sp ecified in the background Minkowski spacetime) 
to (2.1) preserves the Kerr-Schild form, but with trans- 
formed values for H and (see p9|) 
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Now let x a and x'P be the coordinates in the lab frame 
X and the comoving (with the BH) frame X' . For the 
Eddington-Finkelstein system boosted in the xy-plane 
(t>3 = 0), in the lab frame X we have 
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r =7 (Vix + V2y) +x + y +z , 
£ t = t[1 - l{v\x + v 2 y)/r], 
4 = x/r - vxlt, 



(2.6) 
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II. BASIC EQUATIONS 

A. Kerr-Schild form of the Schwarzschild spacetime 

The ingoing Kerr-Schild form of the Kerr metric is 
given by 



ds 2 = (r)^ + 2Rl il i v )dx tl d,x v 



(see |9|,^1(), where /i, v run from to 3, "(]^lv — 
diag(— 1 , 1, 1, 1) is the Minkowski metric in Cartesian co- 
ordinates, H is a scalar function. The vector is null 
both with respect to rj^ v and g^ u , 



rTlJv = .9 M %4 = 0, 



(2.2) 



and we have £ 2 — tli. The spacelike hypersurfaces ex- 
tend smoothly through the horizon, and gradients near 



(2- 1 ) Eq. @ 



where 7 = l/Vl — v 2 and v 2 = v\ + v 2 . x and y are 
defined as x = x — v\t and y = y — v^t. Under a boost 
the metric becomes explicitly time dependent. Since the 
boost of the Schwarzschild solution merely "tilts the time 
axis" , we can consider all the boosted 3 + 1 properties 
at an instant t = Q, in the frame which sees the hole 
moving. Subsequent time t simply offset the solution by 
an a mou nt vt. With Eqs. (2.6) a and Pi are defined via 



B. The scalar field equation 

The field equation for a massless scalar field is 
□ = 0, 
which can be expanded as 



(2.7) 



2 



<r<i> w « = 



-9 



(2.8) 



We can decompose this second-order equation into two 
first-order-in-time equations by defining the auxiliary 
variable it = $ t — f3 i ^ y i, where $ = rip, obtaining, 



$ t = /3 l $,,:+7r, 

7T* = fc + jr. 



(2.9) 



Here 



;r = L ( (3) y$ . . 
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C = ,g M,y (3(^ + V^)(4 + K) - ^ - V^K) 
-?W(^ + ^), 



£ (a 2 (3i t - 2 7 ) - 2 7 - 3/3'V* + 1 
V i 

+ -(l-a 2 ), 
r 

and Vfj, = 7(— is the 4- vector of the velocity v in 
Cartesian coordinates. 

In the comoving (non-boosted) frame, the scalar field 
equation (^]) can be cast into a ID radial equation 
which, for spherical waves, reduces to 



(1 + 2H)d 2 t , + 4Hd t ,d r , + (1 - 2H)d' 2 , - -Hdv 

r 



~Hd r , ~ ^ ) «I> = 0. 



(2.11) 



in-time equations 



Eq. (2.11) can again be decomposed into two first-order- 

(2.12) 
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/3 r = 



r'(l + 2H) r' 2 (l + 2H) 
2H 



(2.13) 
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The cha racteristics (which can b e der ived from metric 
Eq. (|2.l|)) of the scalar field Eq. ( 2.12 ) have speeds —1 
and \r' -2M)/{r' + 2M). 



In the follow ing s ections we solve Eq. (2.9) with a 3D 
code and Eq. ( |2.12| ) with a ID code and compare the 
results. Since the ID code can be used with almost ar- 
bitrary resolution, we can effectively compare our 3D re- 
sults with an "exact" numerical solution. 
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FIG. 1. The initial data in the lab frame X and in the 
comoving frame X' . The bold solid line segments represent 
the wave packets (initial data) in different frames, the dashed 
lines indicate the propagations of the wave packets in space- 
time. The initial data in X (at t = 0) can be obtained by the 
future and past evolution of the scalar field with the initial 
data in X' (at t' = 0). 



III. INITIAL DATA, NUMERICAL ALGORITHMS 
AND BOUNDARY CONDITIONS 

A. Initial value 

As initial data for the scalar field in the comoving 
frame of the BH, we choose a spherical Gaussian of width 
a centered at radius r' 



$(0,x',y',z') = cxp - 



(3.1) 



where = r' + 2Mln(r'/2M — 1) is the tortoise coor- 
dinate (compare (33|). For all calculations in this p ape r 



we choose r' = 10 M and a = 1M. According to (3.1) 



$ vanishes on the event horizon. We also assume time 
symmetry at t = 



dM0,x',y',z')=0. (3.2) 



so that 



n(0,x',y',z') = -P' 8.^(0, x',y',z'). (3.3) 

As the wave packet evolves in time, it splits, with one 
part of it propagating outwards towards null infinity and 
the other propagating inwards towards the horizon. Near 
the horizon, the wave undergoes partial transmission and 
reflection. We calculate the waveform of the scattered 
scalar wave as observed at some fixed distance from the 
BH and compare the 3D results with the "exact" ID 
waveform. 

For the boosted cases, the initial data is different from 
the non-boosted case due to the tilt of axes (see Fig. [l]). 
The frame X'(t', x') (comoving with the BH) moves along 
x-axis with a boost velocity relative to the lab frame 
X(t, x). The initial data in X can be derived by evolving 
the scalar field in X' , followed by a Lorentz transforma- 
tion. (Refer to Appendix |b] for details.) 
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recipe 


advection term 


advection term 


inner 




stable? 






outside BH 


inside BH 


boundary 


v = 


v = 0.2 


v = 0.5 


AI 


C 


C 


E 


Yes 


Yes 


No 


All 


C 


U 


E 


Yes 


Yes 


Yes 


BI 


u 


u 


P 


Yes 


No 


No 


BII 


u 


u 


E 


Yes 


Yes 


Yes 



TABLE I. Summary of recipes AI, All, BI, and BII, and their stability properties for the three different boosts. The boldface 
letters denote C: centered scheme; U: upwind scheme ; E: third-order extrapolation condition; P: copying the time derivative of 
every field at the boundary from the interpolated value just outside the excision zone along the normal direction. We analyzed 
the stability of these recipes for static black holes (v = 0) and boosts with speed v = 0.2 and v = 0.5. Recipe BI, which is based 
on the implementation of |29) passes the non-boosted case but fails in the boosted cases. Recipe AI passes the non-boosted 
case and the v = 0.2 boosted case but fails in the v = 0.5 case. Recipe All and BII pass all tests that we have performed. 



B. Differencing scheme 




♦ boundary points 

centered scheme 

1 | upwind scheme 



FIG. 2. Schematic diagram of shift advection stencils in 
the computational domain and around an excision zone. The 
dashed circles inside the excision zone are the grid poin ts t o 
which the extrapolated values are assigned by using Eq. (3.7). 



We are looking for simple and robust recipes and fo- 
cus, as in p9| , on the shift terms in the field equation, 
which we regard as a generic "advection" terms (terms 
that look like (3 l di). We apply both centered and upwind 
schemes to these advection terms to test their stability. 



All other derivative terms are evaluated using centered 
differencing. 

Consider the gridpoint with array indices i, j, and k 
in the x, y, and z directions, respectively. With centered 
differencing, the first derivative of <f> with respect to x at 
this gridpoint is 



1 



2A; 



[&i+l,j,k — $i-l,j,k]- 



(3.4) 



For the upwind scheme the second-order accurate first 
derivative along x-direction is 



~2Aa 



[$i+2i/i ,j,k 



-4$ 



where v\ is defined as 



for 
for 



P 1 >o, 
/3 1 < 0. 



(3.5) 



(3.6) 



We use analogous differencing in the y- and z-directions 
(see Fig. |). 

Different schemes can be used inside/outside the event 
horizon in each recipe to further isolate where and how 
any instability arises. In this paper, we consider 4 differ- 
ent recipes which we call AI, All, BI, and BII (see Table 
Recipes AI, All, and BII all u se the extrapolation in- 
ner boundary condition (see Sec. [II D| below). In recipe 
AI, the centered scheme is used for finite differencing on 
the shift advection term everywhere. This scheme is sta- 
ble according to a van-Neumann stability analysis since 
/?' is never greater than unity in the Kerr-Schild metric 
(refer to Appendix [X] for a stability analysis). In All, the 
upwind scheme is used inside the BH and the centered 
scheme is used outside the BH for the shift term. In BII, 
upwind differencing is used everywhere. 

Recipe BI, which is based on the implementation of 
p9[ , differs from the other schemes in its excision bound- 
ary condition. Alcubierre and Briigmann p9[ use a cu- 
bic excision zone and copy the time derivative of every 
field at a gridpoint just inside the excision zone from the 
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neighboring gridpoint just outside. Our excision zone is 
spherical with radius r' — M (in the comoving frame) . To 
generalize the scheme of ]29| ] we interpolate neighboring 
gridpoints to the normal on the surface of the excised re- 
gion. More specifically, for a gridpoint (i,j, k) just inside 
the excised region we take its nearest neighbors along the 
coordinate axes away from the center of the black hole, 
say (i + k), + 1, k) and k + 1). These three 
points define a plane, and we interpolate to the intersec- 
tion of this plane with the normal on the surface of the 
excised region. If one of these three neighbors is still in 
the excised region, we project the normal into the plane 
spanned by the remaining two coordinate axes, and do 
the interpolation there; if two of the neighbors are inside 
the excised region we directly copy the remaining third 
point. If all three points are inside the excised region we 
repeat the procedure with the three points + j + fc), 
(i + l,j,k + l) and (i, j + 1, k + 1); if that is also not suc- 
cessful we finally copy the point (i + l,j + l,k+l). This 
particular algorithm is only one of many possibilities. We 
have experimented with a number of other schemes, and 
have found that the stability properties do not depend 
on the details of this implementation. Like BII, BI also 
uses the upwind scheme to compute the shift advection 
terms everywhere. 



C. Outer boundary conditions 

At the outer boundary, we impose an outgoing wave 
boundary condition. In this approximation, we assume 
that the functions are of the form $ = /(At — R) (since 
$ = r<j>), where A = (R - 2M)/(R + 2M) is the outgo- 
ing characteristic speed. The values f(t + At, R) of the 
grid points at the outer boundary are updated by using 
the value f(t, R — AAi) with second-order interpolation. 
In the boosted cases, the characteristics A at the outer 
boundary are derived from the relativistic addition of the 
characteristic in the comoving frame and the boost ve- 
locity. 

This boundary condition provides a stable outer 
boundary provided that the outer boundary is placed at 
a sufficiently large distance. With the outer boundary 
at finite radii, as in our cases, some reflected waves are 
created at the boundary. (We find that reflection waves 
with larger amplitude will appear if the true characteris- 
tic speed A in / is replaced by 1) The amplitude of the 
reflected wave decreases as resolution is increased. In 
boosted cases, the effect of the inaccurate outer bound- 
ary condition cannot be ignored as the BH approaches 
the outer boundary. For our tests, however, the reflected 
waves from the boundary are small perturbations of the 
field and do not affect the stability of different recipes, 
only the accuracy. 



D. Inner boundary conditions 

At the inner boundary we use a singularity excising 
method (see, e.g. @^6|ll|j2|,||-|6)). In general dy- 
namic spacetimes, the location of the apparent horizon is 
not known a priori, and must be computed at each time 
step with an "apparent horizon finder" (e.g. [^7|-^9|). 
This is not the case for the present static Schwarzschild 
background since the location of the apparent (and event) 
horizon is known at all times pQ] . On the other hand, it is 
not necessary to know the exact location of the apparent 
horizon to implement the singularity-excising method, 
provided the excision zone is sufficiently small that one 
can be confident that it lies entirely inside the horizon. 
We choose to excise the region with radius r' — 1M from 
the center of BH where r' is the radial distance measured 
in the comoving frame X'. The simplest inner boundary 
condition obtains, by extrapolation, the field variables at 
gridpoints just inside the masked region from gridpoints 
just outside this region. These values can then be used in 
a centered evolution scheme to update the field variables 
at gridpoints just outside the excision zone. We have im- 
plemented such a boundary condition using a third-order 
extrapolation scheme 



/,-, = 4/ j -6/ j+v +4/ j 



+2i/ 



/ 



(3.7) 



in all recipes except BI. Here j is a gridpoint just outside 
the excised region, and, with v either +1 or —1, j — 
v is just inside (see Fig. ||). We apply this algorithm 
along whichever axis the second derivative is taken, for 
example along the x-axis for a second derivative with 
respect to x. With this extrapolation, the second spatial 
derivatives using centered differencing are second-order 
accurate. This is equivalent to using a second-order one- 
sided differencing scheme. These boundary conditions 
do not violate causality since no information is extracted 
from within the excision zone. This prescription is simple 
to implement and does not require special assumptions 
on the behavior of the variables in the proximity of the 
excision zone. A similar implementation has proven to 
be stable for wave propagation in 2 + 1 dimensions on a 
flat spacetime |?5) . 

In recipe BI, we generalize the implementation of pijj ] 
and copy the time derivative of every field variable at 
the boundary from interpolated values just outside the 
excised region as described in Sec. [II B . As we will show 
in Sec. IV, this inner boundary condition will produce 



stable results in the non-boosted BH case, but results in 
an instability when the BH moves and grid points emerg- 
ing from the excision zone must be assigned by extrapo- 
lation. 



E. Extrapolation to newly emerging grid points 



We use a second-order extrapolation scheme to set val- 
ues at gridpoints that are newly emerging from the exci- 
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FIG. 3. Convergence tests for our ID evolution calculations, using recipe AI (left panel) and BII (right panel). To demonstrate 
second order convergence, we plot the appropriately rescaled finite difference error ErrAr = $Ar/2 — $Ar at the check point 
r — 14M for various resolutions Ar. 



sion zone as the BH travels through the numerical grid. If 
available, we carry out the extrapolation along the z-axis, 
because the BH is boosted in the xy-plane, and the ex- 
trapolated values along the z-direction are least affected 
by numerical errors from previous extrapolations in the 
evolution. If the points needed for this extrapolation are 
themselves excised, we instead average between extrap- 
olations along the x and y directions or use only one of 
the two if the other one is not available (because the 
necessary points are excised). If neither one is available 
we continue through the following hierarchy of preferred 
extrapolation directions: z, then x and y, then xz and 
yz, finally the xy- and xyz-directions. We have experi- 
mented with other extrapolation prescriptions and have 
found that the stability properties of the code are fairly 
insensitive to the details of this scheme. 



IV. NUMERICAL RESULTS 

We numerica lly solve both the ID and 3D scalar 
wave equations (2.9) Eq. ( 2.12 ) using an iterative Crank- 
Nicholson scheme with two corrector steps ffl] . Since the 
ID code can be run with an essentially arbitrary num- 
ber of radial grid points and hence essentially arbitrary 
accuracy, we use this solution as the "exact" solution for 
comparisons. 

We compare the results from the 3D simulations using 
the four recipes (see Table. ||) for both the non-boosted 
and boosted cases. There are two sub cases of the boosted 
case: a "slow" boost speed v = 0.2 and a "fast" boost 
speed v — 0.5. The results for both boosted cases are 
summarized in Table. |, but we will discuss in detail only 
the v = 0.5 boost case. 

In each case, we compare the waveforms with the ID 
result for three uniform grid resolutions A = 0.5M, A = 
0.25M,and A = 0.125M. We choose a time step At 
A I in all cases, and assume equatorial symmetry across 
the z = plane. 



In the non-boosted case, the computational domain 
is 32Af x 32M x 16M (-16M < x < 16A/-16M < 
y < 16M,0 < z < 16M). In the v = 0.5 case, the 
domain is extended to 48M x 48M x 24M (-20M < x < 
28M-24A/ < y < 24M,0 < z < 2AM). A larger spatial 
domain is needed in the high speed case to observe the 
BH's motion before it moves out of the computational 
domain. In each 3D case the code with A = 0.I25M 
resolution is run only to t = 10 M to check its stability 
and convergence; the lower resolution non-boosted cases 
are run until t = 100M, and the lower resolution boosted 
cases are run until t = 48 M. 



A. ID result 

We verify the second-order accuracy of our ID codes 
in Fig. ||. We compare the amplitude difference between 
the results observed from the check point r = 14M by 
continual doubling of the grid resolution. We present re- 
sults for recipe AI in the left panel of Fig. |^ and for BII 
in the right panel. Comparing the scales of ErrAr in the 
two panels shows that, while both scheme arc second or- 
der accurate, the centered scheme of recipe AI converges 
much faster than the more diffusive upwind scheme of 
recipe BII. We find very similar results in the 3D codes 
below. 



B. 3D results in co-moving coordinates 

We first study the four recipes AI, All, BI and BII for 
the non-boosted (comoving) case. All four recipes give 
stable results and converge to the ID result, as shown in 
Fig. § The stability of recipe BI is consistent with the 
conclusion of in which the recipe gives a stable and 
accurate solution in a non-boosted distorted BH. 

In Fig. ^ the first peak and the first minimum in the 
waveform curves are produced by the part of the wave 
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FIG. 4. The scalar field <E> as a function of t at radius r — UM in the non-boosted 3D case. In each panel, the solid curve 
is the ID solution, while the dotted, dashed and dot-dashed line are 3D solutions computed for grid widths A = 0.125M, 
A = 0.25M, and A = 0.5M, respectively. We also include convergence tests for each recipe in the insets, where the dotted 
line is 16 x Erro.i25M, the dashed line is 4 x Erro.25Af, and the dot-dashed line is Erro.25M. Here the finite difference error is 
computed from ErrAa = $ai — $id- 



packet which moves outwards and leaves the computa- 
tional grid. The minimum observed at t « 10M is a nu- 
merical artifact and converges away for increasing grid 
resolution. We demonstrate the second order conver- 
gence of the algorithms in the insets in Fig. ||. For the 
upwind schemes of recipes BI and BII the convergence 
is much slower than for the centered schemes of recipes 
AI and All, as expected from the ID results (see Fig. y) . 
There are several bumps at t « 10M, t w 30M, and 
t 50 M. These are the reflection waves coming from the 
outer boundaries due to the approximate outer boundary 
condition. We find that these bumps can always be di- 
minished significantly by extending the outer boundary 
to a larger radius and/or increasing the resolution. 

The waveforms in Fig. ^ also show that the upwind 
scheme (recipes BI & BII) efficiently eliminates high fre- 
quency noise. This property is very helpful in keeping the 



code stable, but the solution with low resolution can be 
very inaccurate. Therefore, higher resolution and, thus, 
greater computing resources are needed (compared to the 
centered scheme) to achieve a desired accuracy. 

C. 3D results in boosted frames 

We now describe the boosted case with the boost speed 
v = 0.5. In this case the BH, initially located in the 
origin of the lab frame, is boosted along the cc-axis with 
the boost velocity v = (0.5,0,0). We stop the evolution 
when the separation between the center of BH and the 
boundary is less than AM to avoid the BH being too close 
to the boundary. Since the origin of the lab frame is set 
such that the positive boundary of x-axis is at x — 28M, 
the total evolution time is 24M /0.5 = 48M. 



7 



All 




-0.2 

25 50 25 50 25 50 

t/M t/M t/M 



FIG. 5. The scalar field $ as a function of time for a boost with v — 0.5 for checkpoints at (14M,0,0), (0,14Af,0), and 
( — 14M,0,0) (from left to right). The inset in each panel is the blowup of the waveform for early times. The first and second 
rows are the results for recipes All and BII respectively. In each panel, the solid curve is the ID solution, while the dotted, 
dashed and dot-dashed line are 3D solutions for grid resolutions A = 0.125M, A = 0.25M, and A = 0.5M. In the third 
column, for checkpoint (— 14M,0,0), we only show the two coarser resolutions. 



Fig. H shows the waveforms of the scalar field. The 
boosted ID results are derived from the non-boosted (co- 
moving) ID result with a Lorentz transformation (see 
Appendix |b|). We observe the scalar waveforms at 
three check points, located at (14M,0,0), (0,14M,0), and 
(— 14M,0,0). In the first column of Fig. |g, there are two 
peaks in each waveform observed at (14M,Q,0). The first 
peak is the incoming wave packet, and the second peak 
is the outgoing wave packet. This is more obvious in the 
2D snapshots in Fig. || (see also Fig. |l|) . The waveforms 
are masked at t — 26. 3M — 29. 7M in the first column 
of Fig. H because the excision zone passes over the check 
point during that time. The waveforms in the first col- 
umn have larger amplitude changes from low resolution 
to high resolution. The situation is especially obvious for 
recipe BII (using the upwind scheme) in which the peak 



of 4> varies from 0.32 at low resolution to 0.55 at high 
resolution. We restricted the high resolution run to a 
computational domain (-8M < x < 24M,-16Af < y < 
16M,0 < z < 16M) in order to save computational re- 
sources. Accordingly, no high resolution results are avail- 
able for the check point (— 14M, 0, 0). Fig. || shows that 
the 3D results of All and BII still converge to the ID so- 
lution although the rate of convergence for both recipes 
is somewhat worse than in the non-boosted case. 

Recipe BI produces an instability for boosted BHs. We 
believe that this instability arises as a result of combin- 
ing the simple inner boundary condition of recipe BI with 
the treatment of gridpoints that newly emerge from in- 
side the excision zone. The error expands outward as the 
BH moves. Physically, nothing can emerge from the BH, 
but numerically noise can propagate outward through the 
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FIG. 6. The snapshots of the scalar field as evolved by recipe BII in the boosted v = 0.5 case (z = 0.125M). 
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event horizon. Since the third-order extrapolation (3.7) 
condition is equivalent to using a one-sided scheme at 
the excision boundary, it propagates any numerical er- 
rors inwards. We believe that this extrapolation produces 
the i mp roved behavior seen in recipe BII. Furthermore, 
Eq. (3/7) is viable in all cases (non-boosted and boosted) 
as long as the excision zone is inside the event horizon. 

AI is stable in the v = 0.2 case but produces instabil- 
ity in v — 0.5 case. We believe that this instability is 
again caused by the treatment of the grid points emerg- 
ing from inside the excision zone, where numerical er- 
rors are allowed to propagate outwards by the centered 
scheme. The instability can be avoided by using the up- 
wind scheme, which propagate the numerical errors in- 
wards (e.g., recipe All). 

Recipe BII is stable in all the cases. Fig. || shows the 
stable evolution of the scalar field with recipe BII. The 
excision zone can be seen in the snapshots at t = 12M, 
t = 16M and t — 24M where the residual wave ampli- 
tude is seen to hover about the excision zones centered at 
(x,y)=(6M,0), (8M,0), and (12M,0) before disappearing 
entirely at later time inside the masked region. Recipe 
All is also stable in all the cases, although the numerical 
errors in All are bigger than in BII. We have also con- 
firmed that recipes All and BII are stable for boosts in 
off-axis directions by studying a boost of v = 0.5 in the 
v = (0.4,0.3,0) direction. 



grid points emerging from inside the excision zone. We 
find that stability requires a higher-order extrapolation 
condition to update the values of the grid points at the 
boundary of the excision zone when the BH moves. 

Our study may be useful for the numerical evolution 
of binary BHs. It suggests that a simpler scheme (com- 
pared with the causal differencing schemes) might be im- 
plemented to handle evolution and excision for a moving 
BH. Our scalar field model problem provides necessary 
criteria for stability of a numerical recipe. We are now 
gearing up to apply these simple but stable numerical 
recipes to a full dynamical code containing black holes. 
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APPENDIX A: STABILITY ANALYSIS 



V. CONCLUSION 

We have compared four different numerical recipes for 
numerically evolving a scalar field in a Kerr-Schild BH 
spacetime. We have integrated the equations both in 
comoving and lab frames. We study the stability and 
convergence with the goal of identifying simple promis- 
ing techniques for more general evolution problems. We 
used singularity excision and have experimented with the 
boundary conditions for the excision zone. 

Our results show that, in general, an upwind scheme is 
superior to a centered scheme for maintaining stability. 
This is consistent with the results of (^9). We find that 
higher resolution is needed for an upwind scheme than 
for a centered one to achieve a desired accuracy due to 
the diffusive character of an upwind scheme. 

As summarized in Table |, two of the four recipes are 
stable in all testbed scenarios: recipe All, in which the 
upwind scheme is used on the shift (advection) term only 
inside the event horizon, and recipe BII, in which the 
upwind scheme is used on the shift term everywhere, as 
suggested by Both All and BII use an extrapola- 

tion boundary condition for treating the excision zone. 
However, recipe BI, which follows the implementation of 
p9| and simply copies time derivatives as opposed to ex- 
trapolating at the boundary of the excision zone, fails in 
the boosted case. The instability is probably caused by 
the inner boundary condition and its treatment of the 



To analyze the ID radial field equation ( 2.11 ) for sta- 
bility, we simplify it by reta ining only the highest order 
derivative term in / in Eq. (2.12), i.e., 



$,t = /3 r $r+7T, 



(Al) 



where a = 1/(1 + 2H) 2 



1. The Centered Scheme 



We use an iterated Crank- Nicholson method with one 



predictor step and two corrector steps to evolve Eq. ( Al) 
(see ) . To apply a von Neumann stability analysis we 
assume 



<T>" 

j 



(A2) 



and adopt the following finite differencing for the predic- 
tor step ( "forward time centered space" ) 



7T — > 7T 



1 



n+1 



10 



2Aa 



2Ax 



2 sin/cAx 

Ax 
i sin fcAx 



(A3) 



j > 



(Ax) 



(Ax) 2 



(cos fcAx- 1)$™. 



Here n labels the time level and j labels the spatial grid 
point. Substituting (E3) into (Al) yields 



u 



71+1 



At 



Au«, 



(A4) 



where 



A = 



J 



2i/3 r P 1 
2aQ 2i^ r P 

sin fcAx „ 



(A5) 



cos kAx 



2Ax ' ^ (Ax) 2 ■ 
To find the eigenvalues A's for the matrix A, we require 

(A6) 



which gives 
or 

where 
Z 



det(A - AI) = 0, 

(X-2i(3 r P) 2 



2aQ, 



A = 2i 



At' 



(A7) 



(A8) 



At 
Ax 

At sin(fcAx/2 



sin(fcAx/2)(/3 r cos(fcAx/2) ± y/a) 
[2Hcos(kAx/2) ± 1] . 



Ax (1 + 2F) 



(A9) 



Here \Z\ < 1 for the value of the shift (/3 r < 1) and a if 
At/ Ax < 1. 

For the convenien ce of analysis, we would like to decou- 
ple the equation set ( A.4) into two independent equations. 
By using a suitable coordinate transformation R we can 
diagonalize the matrix A 



A' = RAR- 1 = 2i—l, 
At ' 



(A10) 



so that Eq. (A4) decouples into two independent equa- 
tions for the components of the vector u' = Ru, 



u 



ln+l 



u™ = 2iZu™ 



where Z takes two different values (i.e., different choices 
of sign in Eq. (|A9|)) for the two different components of 
u'. 

Now we can use the procedure described in My to de- 
rive the spectral radius (amplification factor) £. The 
first iteration of the iterated Crank-Nicholson method 
start s by calculating an intermediate variable (^u using 
Eq. JA11|) 



(l)~n+l 
3 



(A12) 



where we have dropped the prime for convenience. Av- 
eraging to an intermediate time level n + 1/2 yields 



(l) fi n+l/2 = 



^( {1) u] +1 + u]) = (1 + iZ)u]. (A13) 



This value is now used in the first corrector step, in which 
Wu" +1 ' 2 is used on the right hand side of Eq. (All) 



( 2 )u™ +1 - u? = 2iZ«u" +1/2 = 



= 2iZ(l + iZ)vq. (A14) 

(2) fi n+l/2 = 1 ((2)^+1 + ^ 

= (l + iZ-Z 2 )u]. (A15) 



The second (and final) correc tor st ep now uses ( 2 )u n+1 / 2 
on the right-hand side of Eq. ( All ) 



2 ^(2) fi n+l/2 

2iZ(l + iZ - Z 2 )W*. 



Inserting ( |A2| ) we finally find 



2iZ -2Z 2 - 2iZ 3 . 



(A16) 



(A17) 



Since \Z\ < 1 provided At/ Ax < 1, we find |£| < 1 so 
that the centered differencing scheme is stable. 



2. The Upwind Scheme 

Here we apply a first-order upwi nd s cheme only to the 
shift term. The finite differencing ( |A3| ) won't change ex- 
cept for the terms <!>,, and n r , where 



$ r 



3 ■ = S3> : 



A:> 



Ax 



St: 7 ! 



(A18) 



where S = 



- 1 



Ax 



. Substituting into (|Al[) now yields 



u n+1 - u" 
Zl _L 

At 



= Buy, 



(A19) 



(All) where 
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B 



/ l3 r S 1 
i 2aQ (3 r S 



(A20) 



With the same argument in the last section Eq. ( |A19| ) 
can be rewritten as 



,71+1 



< = 2Wu™ , 



(A21) 



dx^ 



A", 



<9$ 

dx' u 



(B2) 



We obtain it from its definition ir = <& t — (3 l $>.i- By using 
these formula we can derive the initial data (t — 0) and 
the ID waveform (t > 0) in the all boosted cases. 



where 



w = iz - y, 

Y = ^/3 r sin 2 (fcA.T/2) 
Ax 

Z = rhs of Eq.(A9). 



At 2H sin 2 (k Ax/2) 
Ax (1 + 2H) : 



(A22) 



Following the same procedure as in last section one finds 
the spectral radius 



£ = 1 + 2W + 2W 2 + 2W 3 



(A23) 



which again satisfies |£| < 1 provided At/ Ax < 1. The 
first-order one-sided scheme is therefore stable as well. 
For our second-order one -sid ed scheme, we must change 
the finite differencing in (A3) to 



i+2 



4$ 



3+1 



2Ax 



'i+2 



■ 47T 



3+1 



3ttJ 



2Aa 



(A24) 



Analytic analysis becomes complicated in this case, so we 
rely on empirical results, which again exhibit stability. 



APPENDIX B: LORENTZ TRANSFORMATION 
OF A SCALAR WAVE 

In order to produce the initial data for the boosted 
cases as well as the "exact" boosted waveforms to com- 
pare with the 3D numerical results, we need to know 
the transformation of the scalar wave solution in the ID 
static BH background to the solution as viewed in a frame 
in which the BH is boosted. 

Assume a boost velocity v = (vi,V2,0). The Lorentz 
transformation between the comoving frame X' and the 
lab frame X is x ,/3 = A^ a x a , i.e., 



t' = "ft — JViX 



1V2V, 



x H tKvyx + v 2 y) - jvxt, 



y 

z' 



II - 



7 + 1 
7 + 1 



(vix + v 2 y) - jV2t, 



(Bl) 



where 7=1/ Vl — v 2 and v 2 = v\ + v 2 . Since is a 
scalar, it transforms according to $'(x'' 3 ) = $(x Q ). 

For the initial data, we also need to know ir(t = 0,x) 
using 
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